/* 
 * File:   main.cpp
 * Author: jonk
 *
 * Created on 22 de noviembre de 2011, 18:20
 */

#include <cstdlib>
#include <complex>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>
#include<iomanip>
#include <sstream>

#define CONFINED

using namespace std;

#define pi 3.1415926535898

#include "utils.h"

int main(int argc, char** argv) {

    double import(ifstream *MiArchivo, particula *part,  parametros *param);
    void importBinary (ifstream *filein, particula *part,  parametros *param);
    ///Parte de velocidad del tensor de esfuerzos
    double svxy, svxz, svyz, svzz;
    ///Parte del potencia del tensor de esfuerzos
    double spxy, spxz, spyz, spzz;
    ///Tensor de esfuerzos total
    double sxy,sxz,syz,szz;
    ///Energias cineticas y potencial
    double ekin, epot;
	///Tiempo de la simulacion
	double t;
	///Parametros
	parametros param;
	///Particulas
	particula *part = new particula[1000];
	param.u0=50.0/8.0;

	ofstream fileout("stressTensor.dat",ios::out);

	ifstream MiArchivo ("config.bin",ios::in|ios::binary);

	if(MiArchivo.fail())
	{
		cout<<"No se pudo encontrar el archivo config.dat "<<endl;
		exit(1);
	}
	if (MiArchivo.is_open())
	{
	    int nnn=0;
        while(!MiArchivo.eof())
        //while(nnn==0)
	    {

	        nnn++;
	        svxy=svxz=svyz=svzz=0.0;
	        ekin=0.;
            importBinary(&MiArchivo, part,  &param);
            for(int i=0;i<param.npart;i++) {
                ekin+=part[i].mass*(part[i].vx*part[i].vx+part[i].vy*part[i].vy+part[i].vz*part[i].vz);
                svxy+=part[i].mass*part[i].vx*part[i].vy;
                svxz+=part[i].mass*part[i].vx*part[i].vz;
                svyz+=part[i].mass*part[i].vy*part[i].vz;
                svzz+=part[i].mass*part[i].vz*part[i].vz;
            }

            ///Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
            svxy=svxy/(param.vol);
            svxz=svxz/(param.vol);
            svyz=svyz/(param.vol);
            svzz=svzz/(param.vol);
            ///Calculamos la interaccion del sistema sin confinamiento
            epot=interaccion(part,param,&spxy,&spxz, &spyz, &spzz);

            sxy=spxy+svxy;
            sxz=spxz+svxz;
            syz=spyz+svyz;
            szz=spzz+svzz;

            ekin=0.5*ekin/param.npart;

            //fileout<<setprecision(8)<<"  "<<param.t<<"  "<<epot<<"  "<<ekin<<"  "<<sxy<<"  "<<sxz<<"  "<<syz<<"  "<<szz<<endl;
            fileout<<setprecision(8)<<"  "<<param.t<<"  "<<sxz<<endl;
//            cout<<setprecision(8)<<"Sp  "<<epot<<"  "<<spxy<<"  "<<spxz<<"  "<<spyz<<endl;
//            cout<<setprecision(8)<<"Sv  "<<epot<<"  "<<svxy<<"  "<<svxz<<"  "<<svyz<<endl;

        }


		MiArchivo.close();
	}
	delete []part;
	fileout.close();
         return 0;
}



